Effective models for strong electronic correlations at graphene edges 
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We describe a method for deriving effective low-energy theories of electronic interactions at 
graphene edges. Our method is applicable to general edges of honeycomb lattices (zigzag, chiral, 
and even disordered) as long as localized low-energy states (edge states) are present. The central 
characteristic of the effective theories is a dramatically reduced number of degrees of freedom. As a 
consequence, the solution of the effective theory by exact diagonalization is feasible for reasonably 
large ribbon sizes. The quality of the involved approximations is critically assessed by comparing 
the correlation functions obtained from the effective theory with numerically exact quantum Monte- 
Carlo calculations. We discuss effective theories of two levels: a relatively complicated fermionic 
edge state theory and a further reduced Heisenberg spin model. The latter theory paves the way to 
an efficient description of the magnetic features in long and structurally disordered graphene edges 
beyond the mean-field approximation. 

PACS numbers: 81.05.ue,73.22.Pr,73.43.Nq, 73.20.-r 
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I. INTRODUCTION 

The Coulomb repulsion among electrons can give rise 
to exotic phenomena in solid-state materials if this in- 
teraction is sufficiently strong, subsumed under the term 
strong correlations. In graphene, the honeycomb lattice 
of carbon atoms f 1 -^ the question of the strength of the 
Coulomb interaction is, however, not settled yet. On one 
hand there is a strong structural confinement to two di- 
mensions, usually enhancing electronic correlations. But 
on the other hand, the relativistic dispersion near the 
Dirac points with its vanishing density of states (DOS) 
at the charge neutrality point tends to suppress correla- 
tion effects. The bulk of graphene seems to be on the 
verge of the critical interaction strength beyond which 
interaction-driven phase transitions set in^r— 

At the edges of a honeycomb-lattice nanostructure, 
the situation may be significantly different. While the 
strong electronic confinement to a two-dimensional plane 
persists, there is no constraint on the DOS to remain 
small at an edge. Instead, it happens that peaks in 
the DOS emerge, depending on details of the edge ge- 
ometry. Within such spectral ranges, the Coulomb 
repulsion becomes important and leads to interaction- 
induced phenomena, such as edge magnetism at zigzag 
edges^ Edge magnetism has been extensively studied 
theoretically for clean zigzag edges with a great vari- 
ety of methods ranging from self-consistent field theories 
to quantum Monte Carlo (QMC) simulations (see, e.g., 
Refs. [lfj| - [20l ). But edge magnetism is not restricted to 
zigzag edges. Also in chiral ribbons edge magnetism was 
examined theoret ically i 2 1 ' 22 Recently experimental evi- 
dence for spin-polarized edge states has been reported in 
Ref. HH. But it also turned out that in realistic situations 
the magnetic properties of graphene edges depend on the 
environment (passivation, substrate)^ In particular it 
was found that the edge states hybridize strongly with 



the surface state of an Ir(lll) surface, which effectively 
leads to their destruction. 

Most theoretical studies of edge magnetism have been 
performed numerically in the basis of the carbon p z - 
orbitals arranged in a honeycomb lattice. Such calcula- 
tions are typically restricted to relatively small systems 
with a few ten-thousand lattice sites for, e.g., mean-field 
calculations, or even to below thousand sites for quantum 
Monte-Carlo or ab initio methods. 

Here, we follow a different approach. Before we per- 
form the actual calculations, we dramatically reduce the 
degrees of freedom to only the relevant ones, namely the 
edge states. Those have (nearly) zero energy and a strong 
confinement to the edge (in addition to the confinement 
to the 2D graphene plane) . We derive an effective theory 
for edge states at general graphene edges (zigzag, chiral, 
disordered) . For the remaining states (called bulk states) 
it turns out that they can be dropped from the effective 
theory. For observables such as the spin-spin correlation 
function on the lattice, some corrections from the bulk 
states must be taken into account finally. The crucial 
point is that all non-perturbative effects (such as long- 
range correlations) are contained in the effective theory 
for the edge states. The residual bulk states may be ac- 
counted for within perturbation theory or even in a non- 
interacting approximation. Their only effect is a trivial 
local amplification of the magnetic correlations induced 
by the edge states. 

This approximation, i.e., the reduction to an effec- 
tive theory for the edge states only, has been used 
before ) 14 ' 19 ' 20 but in this paper we critically assess the 
quality of this approximation by comparing the results 
of the effective low-energy theory with numerically ex- 
act quantum Monte-Carlo simulations. We find that the 
quantitative accuracy which can be reached within this 
effective theory is at least as good as the limits set by 
uncertainties in the parameters entering the calculations 



2 



(i.e., the hopping amplitude t and the Hubbard parame- 
ter U). The relevant qualitative behavior, i.e., the long- 
range structure of the spin-spin correlations, is remark- 
ably well reproduced. 

Even more important than the high quality of the ap- 
proximation is the fact that the effective theory sepa- 
rates the relevant degrees of freedom from the irrelevant 
ones, and thereby providing valuable insights into the un- 
derlying physics. The corrections from the bulk states, 
mentioned above, actually tend to screen the non-trivial 
edge state effects so that it is sometimes difficult to ex- 
tract the underlying physics, e.g., from the exact results 
of a quantum Monte-Carlo simulation of the full lattice 
model. For instance, the existence of edge magnetism is 
difficult to extract from the exact correlation functions 
for certain chiral ribbons^ since it is not always clear if a 
non-zero correlation function is a bulk effect that decays 
as a power law in the thermodynamic limit, or if it is an 
edge state effect that gives rise to long-ranged correla- 
tions. The effective theory discussed here does not suffer 
from these issues, as it provides direct access to the part 
of the system that potentially leads to long range corre- 
lations and separates this part from the obfuscating bulk 
of the system. 

This paper is organized as follows. In Section HH we 
describe the effective edge state theory and the differ- 
ent contributions to the spin-spin correlation function on 
the honeycomb lattice. In Section IlIIl the spin-spin cor- 
relation function calculated within the effective theory is 
benchmarked against numerically exact quantum Monte- 
Carlo (QMC) results. An intuitive picture of the basic 
mechanisms important for the magnetic correlations is 
developed in Sec. [IVJ Finally, in Sec. E] we perform 
a further approximation, arriving at an effective Heisen- 
berg model for chiral nanoribbons, and show that the 
magnetic correlations in the more complicated fermionic 
theory for the edge states are reproduced. 

II. EFFECTIVE FERMIONIC THEORY FOR 
LOCALIZED EDGE STATES 

A. Derivation of the effective theory 

We start from the nearest neighbor hopping Hamilto- 
nian of a general honeycomb lattice 

H = ' Hc > (!) 

(i,3),r 

where Ci T annihilates an electron with spin r at site i of 
the lattice and runs over nearest neighbors. Let N 
be the total number of lattice sites. Hq can be rewritten 
in terms of its exact eigenstates d^ T = ^ "0^(*) C L' w hh 
Vv(') the normalized single-particle wave function and /i 
a (collective) index, 

H " = X^4t c W; ( 2 ) 

fJ,,T 



where e M denotes the associated eigenenergies. 

The Hubbard Hamiltonian Hjj = U J^i c !t c it c *l c 4 
may be written in the eigenbasis of Hq 

with 

ri234 = E^«^W^a«^ 4 «- (4) 

i 

The central step in the derivation of the effective the- 
ory is the partitioning of the eigenstates of Hq into edge 
states and bulk states 

Kr} = {e pT } U {b VT }. (5) 

For a single zigzag edge this partition is obvious: states 
with zero energy are edge states and those with finite 
energy are bulk states. For arbitrarily-shaped graphene 
structures, however, this energetic criterion is too sim- 
plistic for a reasonable separation of bulk and edge states. 
Instead, one may distinguish edge from bulk states by the 
maximum of the wave function weight = maxj |?/> M (i)| 2 , 
which is a convenient localization measure. If is of or- 
der TV -1 , the state /j is a bulk state, otherwise it is an 
edge state. For a conventional zigzag edge, for instance, 

is of the order of the inverse length of the edge. For 
finite-size systems the partition into edge- and bulk states 
is not rigorous (we will see examples later where it is not 
clear if a certain cigenstate of Hq should be labeled as a 
bulk state or an edge state). In practice, however, it is 
always possible to obtain a feasible partition. If in doubt, 
one may always include more states in the edge-state set. 
In the worst case, if a bulk state is included erroneously 
in the edge-state set, the Hilbert space of the effective 
theory becomes somewhat larger than necessary - the 
physics however remains unchanged. 

With this separation into edge and bulk states, we may 
express the Hubbard Hamiltonian as 

Hu = U Y r i234el t e 2t e^e 4i 

1234 

+ C/^r ' 1234 ]r e\ T e 2T blb if + . . . , (6) 

1234 t 

where the dots refer to terms containing an odd num- 
ber of bulk state operators or four bulk state operators. 
Here and henceforth, an overbar of a binary index, such 
as edge or spin, denotes inversion, i.e. r and f are oppo- 
site spins. In order to shorten the notation, we use the 
same numeric index symbols for edge and bulk states. 
Whether an index, e.g., in a vertex function Ti234 cor- 
responds to a bulk state or to an edge state can always 
be determined from the corresponding operators. For in- 
stance, in Ti234e\ T e 2T el T ,b iT , the indices 1,2,3 run over 
the edge states while index 4 runs over the bulk states. 
If a vertex function or an energy is written without their 
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corresponding electron operators, the bulk/edge indices 
will be indicated by b/e superscripts (r^H^ in the exam- 
ple above). 

Up to now we have only expressed the Hamiltonian 
in a different basis. The first approximation we make 
in approaching an effective low-energy theory is to ne- 
glect all but the first two terms in Eq. We have also 
performed a more controlled approximation based on a 
Schrieffcr- Wolff transformation to eliminate the leading 
order of the terms omitted in Eq. ^ . The resulting addi- 
tional fcrmionic couplings were extremely small, however, 
so that their effect could not be observed in all geome- 
tries discussed in this work. Thus, we have chosen to 
drop those second-order terms completely and keep only 
the terms given in Eq. ([6]). The second approximation 
consists in replacing the remaining two bulk state opera- 
tors by their average with respect to the non-interacting 
Slater-determinant 

b\ T b 2T « 5 12 @(-e b 1 ). (7) 

Thus, the second term in Eq. ^ involves an effective 
edge state hopping t\ 2 e\ T e 2r with 

ti2 = ^E r i e 233©(-4) = tf£[^wr^wp&(i). (8) 

3 i 

where Pb{i) is the electronic density at lattice site i de- 
rived from all occupied bulk states. As was shown in 
Appendix B of Ref. [IH, pb(i) can be calculated directly 
from the total edge state density p e (i) = J2i IV'iWI 2 
by Pb{i) = — Pe(i)) as long as particle- hole sym- 
metry is present.— The site-independent part of pb may 
be dropped as it only leads to a chemical potential 
term which is usually compensated by the positive back- 
ground charge from the lattice. In the particle-hole 
symmetric case, one may furthermore write ^p e {i) — 

S34 "03^)^4 (*)( e 3r e 4r)oj where (-)o is the average with 
respect to the non-interacting ground state in which ex- 
actly half of the edge states are filled. Therefore, the 
effective hopping may be conveniently absorbed into the 
Hubbard part of the effective Hamiltonian 

H eS = E <^ e j«- V + u E r i234 : e t lt e 2f :: e^e 4i :, 

Air 1234 

(9) 

where : e t lr e 2r := e\ T e 2r - (e\ T e 2r ) a . 

The approximations made above are essentially based 
on the assumption that the electronic correlations within 
the bulk states and also the cross correlations of edge 
and bulk states are negligible. In fact, this statement will 
be slightly relaxed subsequently as it turns out that in- 
cluding the bulk background susceptibility enhances the 
agreement of the correlation functions calculated within 
QMC and the effective theory. Nevertheless, the non- 
trivial long-range physics is completely contained in H c h 
as long as U is below its bulk critical strength beyond 
which the whole graphene system (and not only its edge) 
is in an ordered phase. These background corrections only 



lead to an enhancement in the correlation function, the 
basic structure of the correlation function is still solely 
determined by the edge states. 

B. Correlation functions on the lattice 

By now we have derived an effective low-energy theory 
for edge states. For calculating correlation functions on 
the original lattice (this is needed, e.g., for benchmark- 
ing against QMC calculations), we need to translate the 
effective-theory correlation functions back to the original 
formulation on the honeycomb lattice. We consider the 
spin-spin correlation function 

<^> = E rr >t c -sW>> ( 10 ) 

tt' 

formulated in terms of the lattice operators Ci T . These 
may be transformed to a new basis d^ T of Hq eigenstatcs 
so that 

1234 

TT 

(11) 

As in the derivation of the effective theory, we partition 
the d Mr basis states into bulk states b vr and edge states 
e pT . Such a bipartition leads to a large number (2 4 ) of dif- 
ferent terms in Eq. ([lip . However, only few of them give 
rise to significant contributions in the actual correlation 
function. We begin by neglecting all couplings between 
edge and bulk states, for the same reasons as discussed in 
the previous section. Thus, the only non-vanishing terms 
are 

(°i°i}e = E rT '^( i )^(*)V'3 0>4(j)(eLe2r4r' e 4r') 1 
1234 

TT 

1234 

TT 

( a i a j)m = Y^)Mj)Um(j)(b\ T b 2T )(e 3 Ar) 

1234 

T 

+ [eo6], (12) 

and originate from the correlations of the interacting 
edge states (e), the non- interacting bulk states (&), and a 
mixed term (m), involving non-local bulk state and edge 
state densities. In (o~fo-j) m the symbol [e <H- b] stands for 
the first term with all b and e operators interchanged. 

The bulk- and mixed terms will turn out to be im- 
portant only for the sake of a quantitative comparison 
to the QMC correlation functions. They are relatively 
short-ranged compared to the edge state contribution 
and do not contain interesting physics. The edge state 
term, however, contains all the non-trivial effects related 
to long-range correlations and edge magnetism. These 
effects are contained in {e\ T e 2T e^ T ,e 4T')> which will be 
evaluated by exact diagonalization in the present paper. 
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The three terms in the correlation function will be 
shown to be in good agreement with the QMC correla- 
tion functions at graphene edges for small U and repro- 
duce the transition from a local antiferromagnetic cor- 
relation to extended ferromagnetism along the ribbon 
edges. However, if U becomes comparable to t, additional 
bulk corrections originating from the residual interaction 
matrix elements between edge and bulk states must be 
taken into account. This residual interaction on top of 
the edge state correlation function acts as an additional 
Zeeman field to which the bulk states respond linearly. 
This leads to a correction to the correlation function (for 
details see Appendix |A")) 

il 

where 



Xij 



is the bulk state spin susceptibility. The sum only con- 
tains bulk state wave functions. The quality of the cor- 
rection can be further increased by taking the interaction 
between the bulk states into account in the susceptibility. 
This may be done within random-phase approximation 
(RPA), i.e., by replacing 



x b -> x fe ' RPA 



X 



1-Ux b ' 



(15) 



such that we denote the correlation function (7( RPA ) j re- 
spectively. The corresponding background correction in 
RPA reads 

ii 



RPA 6.RPA/ _z z\e 



(16) 



III. THE QUALITY OF THE EFFECTIVE 
MODEL APPROXIMATIONS 



A. The quantum Monte Carlo method 

The effective model is tested against numerically ex- 
act results, obtained by projective auxiliary-field deter- 
minant QMC calculations, in which the groundstate spin- 
spin correlations on finite ribbons are calculated as 



lim 



(V*r|e- 



-8H 



-8H 



(^T\e- 2SH \^ T ) 



(17) 



Here, the trial wave function |?/>t) is required to be non- 
orthogonal to the true groundstate wavefunction. We 
take IV't) as the groundstate of the non-interacting sys- 
tem (U = 0). In the actual calculations, the projec- 
tion parameter 9 is chosen sufficiently large as to en- 
sure convergence to the system's groundstate wavefunc- 
tion. Here, we employed a value of up to 9 = 250, 



which was necessary to meet this requirement. Further- 
more, we implemented a third-order symmetric Suzuki- 
Trotter decomposition of the projection operator, with 
an imaginary-time discretization step At = 0.05, as well 
as an SU(2) spin symmetric Hubbard-Stratonovich de- 
coupling of the local Hubbard interaction term. Further 
details of the algorithm may be found in Ref. [26l . 



B. Effective theory vs. QMC in a zigzag ribbon 

In a zigzag ribbon one edge terminates with A sub- 
lattice sites while the other edge terminates on the B 
sublattice. As a consequence the edge states of opposite 
edges hybridize via the hopping Hamiltonian, while states 
localized at the same edge interact only via the Hubbard 
term. It is therefore convenient to formulate the effec- 
tive theory of a zigzag ribbon in the basis of edge state 
operators efc ST , where k £ [^p ] labels the momentum 
along the edge and s = u,l labels the upper and lower 
edge, respectively, r is a spin label. Note that these basis 
states are not eigenstatcs of the hopping Hamiltonian. In 
this basis the effective Hamiltonian reads 



N 



s k,k' ,q 



(18) 



The primed momentum summation means that 4p < 
k + q, k, k' — q, k' < N x is the length of the rib- 
bon, i.e., the number of unit cells along the edge, tk is 
the hybridization between the edge states of given mo- 
mentum k at different edges. It is equal to the smallest 
eigenvalue of the hopping Hamiltonian for k, which is 
obtained numerically. 

The form of Eq. (jT5J) is remarkable and deserves atten- 
tion. The single-particle term is a momentum-conserving 
inter-edge hopping with hopping amplitude ife. This is 
the only term that effectively couples the two edges; for 
all geometries considered in the present work, the inter- 
edge coupling via bulk states turned out to be much 
smaller than the direct coupling via t^. The Hubbard 
term acts on each edge separately. Thus, the effective 
theory retains the bipartite character of the original hon- 
eycomb lattice. From this very form one may already ex- 
pect that there is an intra-edge ferromagnetic correlation 
mediated directly by the Hubbard term and an antifer- 
romagnetic inter-edge coupling in second order perturba- 
tion theory of the order of ti/U. 

We now discuss the transition from local correlations 
between adjacent edges (small U) to strong ferromag- 
netic correlations along the edges, and demonstrate the 
agreement between QMC results and the effective the- 
ory, which is solved by exact diagonalization. Figure [T] 
shows the intra- and inter-edge correlation functions of 
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FIG. 1. (Color online) Spin correlation functions (ctoct|) at 
the edges of a zigzag ribbon (N x = 7, N y = 6) for three values 
of U = 0.02,0.1,0.3. Site in the first index of the correla- 
tion function refers to site at the lower edge (as indicated by 
the circle). The squares represent the exact results quantum 
Monte-Carlo simulations (error bars are smaller than the sym- 
bol size). The curves show the results of the effective theory. 
The lines between the data points are guides to the eye. The 
dashed red curve is the bare correlation function, consisting 
of edge, bulk and mixed terms. The solid green and dotted 
blue lines (mostly on top of each other) show the susceptibil- 
ity correction without and with RPA, respectively. The edge 
sites are labeled as indicated in the inset. 



the exact QMC results (square symbols) for U < 0.3. For 
U < 0.1 the background corrections are rather small but 
for larger U they become more and more significant since 
the bulk states become locally polarized via the edge 
state correlations (for more details see the derivation of 
the background corrections in Appendix fS]) . However, it 
should be emphasized that these background corrections 
are very short-ranged. They do not affect the long-range 
physics but only enhance the correlations stemming from 
the edge states. Solely the edge states, which are well de- 
scribed by the effective low-energy theory H c ff , determine 
the long-range behavior of the correlation function. 
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a zigzag ribbon with length N x — 7 and width N y = 6 
in the low-t/ regime. We always use periodic boundary 
conditions along the edge. For [7 = 0, the intra-edge cor- 
relation function is exactly a Kronecker delta function, 
while there is a somewhat smeared-out antiferromagnetic 
correlation between adjacent sites at different edges (site 
at the lower edge and site 3 at the upper edge in the 
present example). These antiferromagnetic correlations 
decay along the ribbon edge. As U becomes larger, the 
intra-edge ferromagnetic coupling becomes stronger and 
longer-ranged ferromagnetic correlations start to build 
up within each edge. For larger U each edge hosts a 
super-spin S = N e /2, where N e is the number of edge 
states with different momenta. The inter-edge antiferro- 
magnetic couplings, which have a characteristic strength 
of t\/U become weaker as U is increased, but they are 
still sufficiently strong to produce an antiferromagnetic 
alignment of the intra-edge super-spins in the ground 
state. 

Note that the results of the effective model calculations 
with the most relevant background corrections (bulk sus- 
ceptibility with and without RPA) agree very well with 
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0' 2' 4' 6' 8' 10' 0' 

FIG. 2. (Color online) Spin correlation functions (o-gcrj) at 
the edges of a zigzag ribbon (N x = 11, N y — 8) for larger 
U = 0.2, 1.0. The meaning of the symbols and lines is the 
same as in Fig. [1] 

For large U of the order of the nearest-neighbor hop- 
ping amplitude the background corrections are well vis- 
ible (see Fig. [5]). In the regime where U is larger than 
the inter-edge hopping in the effective theory, the edge 
state contribution to the correlation function is saturated 
and thus becomes independent of U. In this regime, 
only the background correction increases, but the edge 
state correlations are saturated. Also, for larger U there 
is a visible difference between the non-interacting back- 
ground correction C^ 1 ' and the RPA background correc- 
tion C( RPA ' . Naturally, if U comes too close to the RPA 
critical U c ~ 2.23^ the background correction is overes- 
timated. This is why it is difficult to reproduce the exact 
QMC correlation function for that large U as this would 
imply the need for high order corrections in perturbation 
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theory. However, for understanding the underlying long 
range physics, which is contained in the edge state sub- 
system, it is not even necessary to reproduce the exact 
lattice correlation function with all corrections sitting on 
top of the actually relevant edge state physics. The study 
of the edge contribution plus all corrections rather obfus- 
cates the simple physical picture that will be discussed 
in the next section. 



C. Effective theory vs. QMC in chiral ribbons 

The derivation of an effective low-energy theory is not 
restricted to zigzag ribbons. Also for chiral ribbons we 
may obtain an effective edge state theory along the gen- 
eral lines described in Sec. Ill Al We characterize the rib- 
bon geometry by the number of unit cells along the rib- 
bon N x , the number of zigzag lines across the ribbon N y , 
the length of one unit cell x, and an additional shift 6 of 
the unit cells. Figure [3] clarifies those definitions. Note 
that the shift 6 changes the chirality of the ribbon in a 
non-standard way (standard would be © = 0). We intro- 
duce 6 in order to be able to tune the localization of an 
edge state along the edge (see below). 

The eigenstates of Hq may be obtained in fc-space 



ij)k a {m,n) = 



1 



l 4>ka(n), 



(19) 



where m runs over different unit cells, n runs through the 
sites within one unit cell, and 4>ka(n) is the transverse 
wave function for the ath eigenstate with momentum k. 

In Fig. 2] we show the band structure of a standard 
(6 = 0) ribbon with \ = 4 and N y = 200. Obviously 
there is a clear partitioning between edge states and bulk 
states in this chiral ribbon. In ribbons with different 
parameters one may find similar partitions, although the 
edge states may not exist for all momenta in the Brillouin 
zone. 



(a) 

X = 4 rrYVY 


(b) ^Y^Yy 
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FIG. 3. (Color online) Definition of the shifted chirality ge- 
ometry via x, 6, and N y . Dashed lines indicate connections 
between sites of neighboring unit cells. Part (a) shows a stan- 
dard chiral nanoribbon. Part (b) shows a nanoribbon with an 
additional shift & — 1 of the unit cells. 

As for the zigzag ribbons we compare the spin corre- 
lation function (ofcj) for i, j right at the edge. We do 
so for two chiral geometries: in one case (geometry 1) 
we use N y — 8 and 6 = (Fig. [S]) and in the other 
case (geometry 2) we use N y = 10 and 6 = 1 (Fig. [5]). 




FIG. 4. (Color online) Edge states and bulk states of a chi- 
ral ribbon with \ — 4 and N y — 200. The inset shows the 
geometry along one edge. 



N x = 5 and x = 4 in both cases. As far as the agreement 
of QMC and the effective theory is concerned, the result 
is very similar to what we have found for zigzag ribbons. 
For small U = 0.2 the agreement is perfect and the back- 
ground correction is very small. For larger U = 1 we find 
a considerable background correction that enhances the 
underlying edge state correlations. 

More important than a perfect reproduction of the 
QMC results by the effective theory is the qualitative 
agreement with respect to the suppression of the fer- 
romagnetic order along the ribbon. In case of geome- 
try 1 the ferromagnetic correlation is already strong for 
U = 0.2 and for U = 1 the edge state ferromagnetism 
is essentially saturated and the correlations are extended 
along the edges. For geometry 2, however, the ferromag- 
netism along the ribbon is strongly suppressed even for 
U = 1. Wc attribute this suppression to the weakened 
ferromagnetic interactions due to 6 > 0. It is most im- 
portant to notice that the effective model agrees with the 
exact QMC solution regarding even this subtle point. 



IV. AN INTUITIVE PICTURE 

The bipartite nature of the honeycomb lattice trans- 
lates to a bipartiteness in the effective edge state the- 
ory. The spectrum of Hq is particle-hole symmetric at 
half-filling, i.e., all eigenstates come in pairs of positive 
and negative energies ±e Q . Let their wave functions be 
(j) a ±(i) £ IL The two linear combinations </> a+ (z)+</> Q _(i) 
and <j) a +(i) — </> Q _(z) live on different sublattices. This 
principle obviously holds also for the edge states. Thus, 
we may choose the edge state basis e asT where s — A, B 
labels the sublattice, on which the wave function lives. 
One example for this sublattice-resolved basis was used 
in Sec. IIII Bl where a was a momentum. Another ex- 
ample will be studied in the next section, where a labels 
the center of a Wannicr state. In this sublattice-resolved 
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FIG. 5. (Color online) Spin correlation function (a^tJj) at 
the edges of a chiral ribbon for two interaction strengths U = 
0.2, 1. The number of unit cells is N x = 5 and the number 
of zigzag lines within each unit cell is N y = 8. The positive 
correlations are intra-edge and the negative correlations are 
inter-edge. The edge sites labeling is indicated in the inset. 
The meaning of the symbols and lines is the same as in Fig. 

hi 



0.02 



0.01 



10 



20 



U = 0.2 




FIG. 6. (Color online) Spin correlation function (cr|(j|) at the 
edges of a chiral ribbon with shift 6 = 1. The number of unit 
cells is N x = 5 and the number of zigzag lines within each 
unit cell is N y = 10. The positive correlations are intra-edge 
and the negative correlations are inter-edge. The edge sites 
labeling is indicated in the inset. The meaning of the symbols 
and lines is the same as in Fig. Q] 



basis, Hq is of the form 

H = ^2 taa'cisrCa'ST, (20) 
saa'r 

i.e., it couples only states on different sublattices. For 
clean zigzag ribbons the edge states on different sublat- 
tices also live on different edges. In contrast, the Hubbard 
Hamiltonian 

H u = U ^2 ri234, s eJ st e 2s te3 S ; e 4sl (21) 

1234, s 

couples only states on the same sublatticc, as is obvious 
from the definition of the effective vertex F in Eq. (0}. 

Based on this bipartitioning of the edge states into 
groups on different sublattices one arrives at the follow- 
ing picture: Hjj essentially gives rise to a ferromagnetic 
coupling of order U within each sublattice. On the other 
hand, Hq gives rise to an antiferromagnetic coupling of 
order \t aa '\ 2 /U between different sublattices. This pic- 
ture becomes especially clear if instead of plane waves 
along the edge, one uses a localized Wannier basis. This 



will be discussed in detail in Sec. |V] For perfect zigzag 
ribbons it turns out that the fc-space basis is superior to 
the localized Wannier basis. But for chiral ribbons and, 
as we expect, also for disordered ribbons, where the edge 
states are naturally more localized along the edge, the 
Wannier basis turns out to open up promising possibili- 
ties for a further simplification of the fermionic edge state 
model to a Heisenberg model. 

For the zigzag edge states, however, we may still 
demonstrate the qualitative validity of the intuitive pic- 
ture on a heuristic level. For this we investigate the ex- 
citation spectrum of the edge state subsystem by exact 
diagonalization of the Hamiltonian of the fermionic ef- 
fective theory. The excitation energies of the bulk states 
are much higher so that they will not interfere in this 
low-energy consideration. In these spectra we will iden- 
tify the traces of the two different couplings (intra- and 
inter- edge). 

Figure [7] shows the lowest four excitation energies of 
a zigzag ribbon with length N x = 11. Such a rib- 
bon has four edge states at each edge with momenta 
k = 871-/11, IOtt/11, 127T/11, 14tt/11. For large U the four 
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u 

FIG. 7. (Color online) Lowest excitation energies of a zigzag 
ribbon (N x = 11, N y = 20) as a function of U. For orientation 
the ground state energy (E = 0) is included in this plot. 
The dots are results from the exact diagonalization of the 
effective edge state Hamiltonian. The lines show the super- 
spin approximation for large U (see text). 

edge states at each edge form 5 = 2 super spins Sa and 
Sb, since the intra-edge coupling mediated by Hjj [Eq. 
(f2"Tj) ] is dominant in this regime. The effective inter-edge 
hopping H [Eq. (f20|) ] couples the super spins antiferro- 
magnetically, which may be expressed as JSa ■ Sb- The 
total spin J = Sa + Sb is zero in the ground state and 
the excitation energies are + 1), with j = 1, . . . , 4. 
In order to show that the inter-edge coupling J is in- 
versely proportional to U, we fit + 1) to the large 
U part of the exact diagonalization results and get very 
good agreement for large U - not only for single curves 
with fixed j, but also for different excited states. For 
small U the ED results deviate from the super-spin ap- 
proximation, which signals the breaking up of the super 
spins. 

At this point it is easy to understand the saturat- 
ing behavior of the edge state theory for large U . For 
small U the antiferromagnetic coupling of edge sites at 
opposite edges is relevant and thus suppresses the long 
range ferromagnetic order along the ribbon. As U is 
increased, the intra-edge ferromagnetic correlations be- 
come stronger. In the limit of large U each edge (of fi- 
nite length) is ferromagnetically ordered and constitutes 
a rigid super-spin. The much weaker antiferromagnetic 
inter-edge coupling then forces the super spins into an 
antiparallcl alignment, but is not strong enough to de- 
stroy the super spins and therewith the ferromagnetic 
ground state correlations that are extended along the 
(finite-sized) ribbon edge. 



V. EFFECTIVE HEISENBERG MODEL FOR 
CHIRAL RIBBONS 

Based on the analysis given above it is tempting to de- 
velop an intuitive picture on the basis of effective spins 



at the (zigzag or chiral) edge, being ferromagnetically 
coupled along the same edge and antiferromagnetically 
coupled between different ribbon edges. We now demon- 
strate that, in order to actually obtain a satisfactory de- 
scription in terms of effective spins, one needs to achieve 
a certain degree of spatial separation of the effective edge 
spins in real space. For zigzag edges the edge states are 
not sufficiently separated for a quantitative agreement. 

A. Wannier edge states 

The effective theory for edge states at perfect zigzag 
ribbons is best formulated in momentum space. It is of 
course possible to transform any momentum space basis 
to a real space basis, i.e., a maximally localized Wannier 
basis (henceforth we will call such a localized basis Wan- 
nier edge states). However for zigzag edges the Wannier 
edge states arc rather dclocalizcd along the edge so that 
neighboring wave functions have a considerable overlap. 

In chiral ribbons, where the edge consists of a scries 
of zigzag segments separated by steps (see Fig. [3]) , it is 
possible to obtain improved Wannier edge states. Each 
zigzag segment hosts one localized edge state and each 
step acts as a barrier between neighboring zigzag seg- 
ments. In order to be able to further tune the Wannier 
edge state separation, we will make use of the shifted 
chirality geometry, in which the barrier height & may be 
changed by shifting the zigzag segments relative to each 
other. Zigzag segments of length x = 4 are especially 
convenient as a starting point, since in the folded Bril- 
louin zone (see Fig. [2J there is an edge state for each 
momentum k along the edge (remember that the edge 
state momenta are usually restricted to some interval, 
e.g. [27r/3,47r/3] for zigzag edges). 

Given the transverse wave functions <fiks (n) of the edge 
states, where n labels the site within one unit cell and s 
labels the edge, the Wannier edge states are 

Vw(m,n) = ^-J2 c ~ ik{X ' m) ^s(n). (22) 

k 

Here, x is an integer labeling the unit cell to which the 
Wannier edge state is localized. <fiks(n) is obtained nu- 
merically from the hopping Hamiltonian Hq in k-space, 
and so is the hybridization amplitude tk between <pks and 
cj)ks- The arbitrary phase of (j>k s {n) is fixed by requiring 
Em n IV'a:s( m > n )\ A to be maximal. Thus, each pair x, s 
corresponds to one spin-degenerate localized state. In 
what follows we will see that the largest part of the Hub- 
bard interaction will force one single electron to such a 
state, thus generating a localized spin. 

A convenient measure for the degree of localization is 
the ratio 

Ri Em,J^,s(m,n)| 4 

Em,n \4>x,a(m,n)\ 2 \4>x+l,s( m > n )\ 2 

of the self-overlap and the overlap of neighboring Wan- 
nier edge states. The numerator of i?i oc corresponds to an 
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effective Hubbard U* , i.e., the energy penaity of a double 
occupation of one Wannier edge state. The denominator 
will later be interpreted as the nearest-neighbor part of 
the ferromagnetic intra-edge coupling. 

For a zigzag ribbon, i?i oc is typically between 4 and 5 
(depending on the length and width of the ribbon). For 
chiral ribbons with variable & the behavior of i?i oc with 
system size is illustrated in Fig. [8l It grows dramatically 
as 6 is increased. Since the Heisenberg model becomes 
better for larger i?i oc it is clear that the Heisenberg ap- 
proximation is better for larger 6, i.e., for larger barriers 
along the edge. 



100 




FIG. 8. (Color online) _Ri oc as a function of the ribbon geom- 
etry N y = 2N X and ©, as indicated. 

The effective edge state Hamiltonian can be expressed 
in the basis of Wannier edge states. The hopping part 
reads 



tf h op = £ t xxl 



^XST^X' ST H.C. , 



with 



t / — — X^fj J k ( x - x ') 

^XX 1 j. j / L K C 



(24) 



(25) 



-ffhop couples Wannier edge states at different edges. In 
contrast, the effective interaction is only non-zero if all 
states participating in one vertex are at the same edge 

Hjj = U F XlX2X3Xi:S e Xis ^.e X2S fe X3S ^e Xi si- (26) 

1234, s 



B. Heisenberg model 

As the Wannier edge states becomes more and more 
separated, the edge states increasingly behave like a se- 
ries of Heisenberg spins. We now derive the effective 
couplings between those edge spins. For this we consider 
pairs of localized Wannier basis states and perform an 
up to second order perturbation theory expansion with 
l/UF xxxx as the small parameter. 



We start with two states e x with x = 1, 2 at the same 
edge s, therefore sitting on the same sublattice. The 
effective hopping between states on the same sublattice 
vanishes and we are left with three terms from Hjj (we 
drop the index s) 

Hit = U ^ ' r ' X xxx e x f e xf e x ± e xl i (27) 

3=1,2 

H^P = UT U22 YSAr^Af^r - e\ T e lf e\ f e 2 r) , (28) 

T 

H' = C/r 1 i22e} t e 2 te{ 4 _e24. + U £ Fv2a3,Ar e ^Af e ^ 

T 

+ U^T 2111 el T e lT e\ f e 1 r + H.c. (29) 

T 

The first term is a Hubbard interaction for the lo- 
calized edge states. It forces the states 1 and 2 to be 
occupied with one electron, respectively. Since Hy is 
the dominant term in the Hamiltonian it is a good ap- 
proximation to restrict the Hilbcrt space to those many- 
body states having one electron per Wannier edge state. 
Within this restricted Hilbert space H { p acts as a ferro- 
magnetic coupling between the spins of the electrons in 
the two Wannier edge states. H' has no non-zero matrix- 
elements in this reduced space. In principle it is straight- 
forward to calculate the second order perturbation theory 
correction due to H', i.e. — i/'(iJ^) _1 iJ'. However, for 
the geometries we are considering here this correction is 
extremely small and we were not able to find any no- 
ticeable effect of these terms on the correlation function. 
In the following we will drop this correction. Thus the 
intra-edge ferromagnetic coupling can be written as 



H 



FM 



x<x' ,s 



J. 



XX' xs 



UT 



XXX' X' ,8 



(30) 

where cr xs is the vector of Pauli matrices corresponding 
to the spin of the electron in Wannier edge state x at 
edge s. 

Next we consider a pair of Wannier edge states e x with 
x = 1,2 at different edges, therefore sitting on differ- 
ent sublattices. In this case there are only two terms 
in the Hamiltonian. The dominant term is again 
so that we again may restrict ourselves to the case of 
one electron per edge state. The second part is the ef- 
fective hopping -f/hop [Eq. (|24p] which has no matrix 
elements in the restricted subspace and must therefore 
be taken into account in second order perturbation the- 
ory via — -ffhop(-fftj) ^^hop- ln this case one obtains the 
standard antifcrromagnetic interaction between the spins 
of the electrons in edge states 1 and 2, i.e., 



■HaFM = £ Jxx' a xs 



with 



TSS 

fJ xx' 



t xx i{y xxxx.s ~t~ ^ x'x' x' x' ,s) 
^UTxxxx s^x'x'x'x 1 - 



(31) 



(32) 



10 




1' 2' 3' 4' 5' 

Wannier state 




1' 2' 3' 4' 5' 

Wannier state 




1' 2' 3' 4' 5' 

Wannier state 



FIG. 9. (Color online) Effective correlation function {afa%) 
of the Wannier edge states of chiral ribbons. The parameters 
are: N x — 5, N y — 56, and & — 8 in part (a); N x = 5, 
N y = 120, and 6 = 8 in part (b); iV* = 5, N y = 10, and 
& = in part (c). The insets show the lattice geometry to- 
gether with the positions of the Wannier states. U = 2 in all 
graphs. The reference site is marked by a red circle. The pos- 
itive correlations correspond to Wannier states at the same 
edge as the reference site. The dots are calculated from the 
fermionic edge state theory. The line connects points calcu- 
lated from the Heisenberg model, which is an approximation 
of the fermionic theory. The line segments connecting the 
Wannier state numbers are only guides to the eye. 



Note that this expression also holds for edges with- 
out translation symmetry (e.g. disordered edges). For 



the chiral ribbons with translation symmetry we have 
UT XXXXS = U* for all x, s and thus the expression for 
the antiferromagnctic coupling reduces to the well-known 

^xx' = ^xx'/U ■ 



C. Fermionic theory vs. Heisenberg theory 

In order to scrutinize the Heisenberg approximation 
we compare the correlation functions of the spins in the 
singly-occupied Wannier edge states, calculated within 
the fermionic edge state theory and in its Heisenberg ap- 
proximation. In both cases the ground state is calculated 
by exactly diagonalizing the effective Hamiltonians. The 
fermionic correlation function is 

tt' 

(33) 

The results of this comparison for three example ge- 
ometries and U — 2 is shown in Fig. [9l In part (a) and 
part (b) the Wannier edge states are well separated due 
to large chiral shifts 6 = 8, which results in a nearly 
perfect agreement between the two edge state theories. 
In part (c) a shift (3 = geometry is considered. In 
this case, there are apparent differences between the cor- 
relation functions calculated from the fermionic theory 
and from the Heisenberg approximation, which can be 
traced to the fact that the neighboring Wannier edge 
states have a relatively large overlap (i?i oc ~ 8.4 in this 
case) . Still the qualitative behavior of the magnetic cor- 
relations along the edges is reproduced. 

An important point to note is that the effective Heisen- 
berg theory is able to reproduce the transition between 
the local inter-edge antiferromagnctic correlation [part 
(a) of Fig. [3] and the extended intra-edge ferromagnetic 
correlation [part (b) of Fig. [9]. The exact solution of 
the fermionic edge state theory via exact diagonalization 
is limited to O(10) Wannier edge states. The effective 
Heisenberg spin model, however, is solvable for several 
thousand spins by means of highly efficient world-line 
QMC simulations, even in the presence of extended ex- 
change interactions^. The crucial observation here is 
that for the effective Heisenberg model no QMC sign 
problem emerges, because the bipartiteness of the honey- 
comb lattice underlying the original Hubbard model de- 
scription translates into a bipartitioning of the effective 
Heisenberg spin model, along with ferromagnetic (anti- 
ferromagnetic) Heisenberg exchange interactions among 
sites that belong to equal (different) sublattices. This 
remarkable feature reminds of a similar commensura- 
bility effect in the RKKY-interactions among magnetic 
adatoms or lattice-defcct-induccd local moments medi- 
ated by the graphene bulk electrons on the honeycomb 
latticed. 
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VI. SUMMARY AND DISCUSSION 

We have shown that the magnetic features of graphcnc 
edges that are due to electron-electron interactions can 
be described solely on the basis of effective theories for 
the edge states. We have given simple and general rules 
for the construction of these effective theories, applicable 
to arbitrary edge geometries. The observables defined 
on the original honeycomb lattice, the most important 
of which is the spin-spin correlation function, can be re- 
constructed from the edge state correlation function pre- 
dicted by the effective theory. This reconstruction in- 
volves a background correction from the bulk states, but 
only in a non-interacting or mean-field approximation, so 
that no elaborate methods are needed for its evaluation. 
In fact, one obtains already a reasonably good estimate 
of the background correction by simply multiplying the 
edge state correlation function from the effective theory 
by a factor of 1.5U. This shows that the background 
correction has essentially the effect of a trivial local am- 
plification of the spin-correlations arising from the edge 
states. 

Based on this observation we argued that in order to 
understand the correlations along graphene edges it is 
advantageous to study the effective theory directly. The 
main reason for this is that the complicated structure 
of the lattice-resolved correlation function tends to ob- 
scure the underlying physics. This can be seen easily 
by comparing Figs. [5] and [9l describing essentially the 
same physics, but in the first plot the original lattice is 
reconstructed, while in the second plot it is not. 

If the edge states can be written in a real space basis 
with well separated wave functions, a further approxima- 
tion of the fermionic theory, namely a Hciscnbcrg theory 
of edge states, is feasible and able to reproduce the effec- 
tive correlation functions. We have given simple formu- 
las by which the Heisenbcrg coupling constants can be 
evaluated directly from the wave functions of the hop- 
ping Hamiltonian. Furthermore, we have shown that the 
fermionic edge state theory (which agrees well with the 
exact solution) agrees even quantitatively with the much 
simpler Hcisenberg theory for the edge states, as long as 
the edge states are well separated in their maximally lo- 
calized Wannier basis. Even for less well separated edge 
states we still find qualitative agreement. 

The effective theories described here enable the theo- 
retical study of graphene systems of realistic sizes. Espe- 
cially with the Hcisenberg model it is possible to study 
thousands of spins, which corresponds to ribbon lengths 
of micrometers. 
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Appendix A: Background correction 

The exact zero temperature spin-spin correlation func- 
tion between two sites i,j on the honeycomb lattice reads 

( CT ?a*> = c(i,j) = E Tr ' E r^M^mMj) 

tt' 1234 

x (^ \d\ T d 2 rdl T ,d iT ,\^ ), (Al) 

where d^r annihilates an electron with spin r in a state 
(i with wave function ^^(i). We assume that d^ T is an 
cigenstate of Hq [Eq. @] with eigenvalue e^. |^o) is the 
ground state, the average with respect to which may be 
written as (following the standard procedure described 
in Ref. M) 

(*0M t lr d 2r4r' rf 4r'|*0> 
°° / ;\n poo 

x ($o\T[H 1 (t 1 )...H 1 (t n )dl T d 2T dl T ,di T ']\$o}c, (A2) 

where H\(t) is the interaction picture operator of all 
terms of H\j [Eq. ([3])] but the ones involving four edge 
state operators. |$o) is the ground state of H . The 
terms with four edge state operators are accounted for 
exactly in the effective edge state theory. Note also that 
only connected diagrams are to be included in this series. 
We may write Hi (t) as 

J?i(t) - L/^'r 12 34d t lt W d 2tW4|Wrf4i(i), (A3) 

1234 

where the primed sum means that terms with four edge 
state operators are excluded. Since in the only terms 
coupling edge and bulk states are contained in Hi, the 
bulk state operators may be contracted separately in Eq. 
(|A2j) . while the edge state averages are left to be calcu- 
lated exactly. 

It turns out to be sufficient to only retain the first or- 
der in U in the perturbation series (|A2j) . From the corre- 
sponding average (§o\T[Hi(t)d\ T d2 T dl T ,d4 T r)\Qo) c we se- 
lect only the terms in which there is one pair of edge state 
creation/annihilation operators in Hi and one pair in the 
group d\ T d2 T d\ T ,d4 :T ' . We drop all other terms since they 
are smaller and the quality of the background correction 
reached by this lowest order correction is already much 
better than the error bars involved in the parameters 
(the Hubbard U used in the literature fluctuates by fac- 
tors up to 4) entering the initial model. Besides it is not 
our aim to actually reconstruct the exact result as good 
as possible. We rather argue that it is in principle possi- 
ble to increase the agreement with the exact correlation 
function in a perturbative way, but for understanding the 
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underlying physics it is not recommended to work in a 
lattice formulation at all. 

I 



Within the approximation described above the first or- 
der correction to the correlation function reads 



h 

(-i) / dte-"l*l(e^ i e 6Tie t lT e 2T )($o|T[6^ i (t)& 8fl (t)4, &4T ,]|$ ) c . (A4) 



1...8 

tt'ti 



Remember that in our notation the operator structure 
determines whether the numeric indices 1 ... 8 run over 
bulk or edge states. In this expression we calculate the 
edge state correlation function {e\ Ti eQ Tl e\ T e2 T ) from the 
effective edge state theory. The remaining four point 
function involving the bulk state operators may be ex- 
pressed in terms of the bulk state spin susceptibility 



8 unocc. 
7 occ. 



£8 - £7 



(A5) 



where the state summation is restricted to the bulk states 
only. A straightforward calculation then gives 



(A6) 



This expression for the background correction gives 
rise to an appealing interpretation. For this we assume 
that i and j are edge sites with a distance of, say, at least 
5 lattice constants. One may understand the measure- 
ment of the spin-spin correlation function between two 
sites i and j as fixing the spin at site i and then mea- 
suring the mean spin at site j. Fixing the spin at i leads 
to a long-range spin polarization coming from the edge 
states, expressed within the effective edge state model 
by (oferf^) 6 . This polarization may be measured directly 



I 

at j by setting i\ = j. But there is also the effect that 
the bulk states, which are themselves not capable of de- 
veloping long-range correlations, become spin-polarized 
near a non- vanishing edge state polarization at i\. This 
is because the bulk states feel an effective Zeeman field 
7^(of ofj 6 at site i\ due to the edge state spin correla- 
tions. This additional spin polarization of the bulk states 
at site j due to the edge state polarization at site i\ is 
described by the bulk spin susceptibility x\ x j- This is the 



first term in Eq. (|A6[) . The second term is the symmetric 
process under the exchange of i and j. 

For the geometries studied in the present paper, x% 
depends separately on i and j. However, it is always true 
that the diagonal contribution Xu — 0.4 ±0.1 is by far 
the largest. Thus, the main effect of taking the correction 
C^'(i,j) into account amounts to a multiplication of the 
edge state correlation function by (1 + 0.4/7). 

Following the argumentation above, higher order (in 
U) corrections may be taken partially into account by re- 
placing Xij by its RPA series. As shown in the main part 
of the paper, doing so increases the agreement with the 
exact correlation function for U below the mean-field crit- 
ical interaction for the antifcrromagnctic instability. For 
U close to Z7 cr it., however, the overestimation of the anti- 
ferromagnetic spin response by RPA naturally leads to a 
considerable overestimation of the total spin-spin corre- 
lation. Nevertheless, the qualitative long-range structure 
of the spin correlations can be determined solely within 
the effective theory for the edge states. 
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